#################################################################################### 
####################################################################################
####################################################################################
# Appendix E
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
library(ggplot2)
####################################################################################
####################################################################################
## Read in survey data (from: analysis_results.R)
dat = readRDS(here("data/dat_unrestricted.rds"))

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Bushenyi","Ibanda","Jinja","Mbale","Mpigi","Pallisa","Mbarara","Shema")
dat[, paste0("c_", cols)] = lapply(dat[, cols], function(x) (x - mean(x))  )
covariates = c(#"sex + age + education_level + Arua + Bushenyi + Ibanda + Jinja + Mbale + Mpigi + Pallisa + Mbarara",
  #"Arua + Bushenyi + Ibanda + Jinja + Mbale + Mpigi + Pallisa + Mbarara",
  "treatment*c_sex + c_age*treatment + c_education_level*treatment + c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment", 
  "c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment")


####################################################################################
####################################################################################
## Social Desirability Bias ####

sdb <- dat %>% 
  group_by(sent_us) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

sdbp = ggplot(sdb[with(sdb,order(-perc)),], aes(x = reorder(sent_us,-perc), y = perc*100)) +
  geom_bar(stat="identity") +
  scale_x_discrete(labels= c("Researchers", "Government", "An NGO", 
                             "Other", "Don't Know")) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(sdb$perc*100)+1.5), labels=function(x) paste0(x,"%")) +
  labs(title = NULL,
       x = NULL,
       y = NULL) + 
  theme_bw() +
  geom_text(data=sdb[with(sdb,order(-perc)),], aes(label=paste0(round(perc*100,1),"%"),
                                                   y=perc*100-3), size=4, family="CM Roman") +
  theme(axis.title = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(hjust = 0.5),
        text=element_text(family="CM Roman")) 

rm(sdbp,sdb)

## Who Sent Us?
for (i in names(dat[,c("sent_us_ngo","sent_us_govt","sent_us_ir","sent_us_dk")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

stargazer(m.sent_us_ngo.c, m.sent_us_ngo.t, m.sent_us_govt.c, m.sent_us_govt.t,  
          m.sent_us_ir.c, m.sent_us_ir.t, m.sent_us_dk.c, m.sent_us_dk.t, 
          
          se = list(c.sent_us_ngo.c[,2], c.sent_us_ngo.t[,2], c.sent_us_govt.c[,2], c.sent_us_govt.t[,2], 
                    c.sent_us_ir.c[,2], c.sent_us_ir.t[,2], c.sent_us_dk.c[,2], c.sent_us_dk.t[,2] ),
          
          p = list(c.sent_us_ngo.c[,4], c.sent_us_ngo.t[,4], c.sent_us_govt.c[,4], c.sent_us_govt.t[,4], 
                   c.sent_us_ir.c[,4], c.sent_us_ir.t[,4], c.sent_us_dk.c[,4], c.sent_us_dk.t[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/social_desirability.tex",
          label = "tab:social_desirability",column.sep.width = "1pt",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of CHP Intervention on Beliefs About Researchers",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("An NGO","Government","Researchers","Don't Know"), column.separate = c(2,2,2,2),
          add.lines = list(c("Covariates", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

rm(list =  grep("^(c|m)\\.", ls(), value = T))

####################################################################################